C     双椭球热源模板
      SUBROUTINE DFLUX(FLUX,SOL,JSTEP,JINC,TIME,NOEL,NPT,COORDS,JLTYP,
     \                 TEMP,PRESS,SNAME)

      INCLUDE 'ABA_PARAM.INC'

C     TIME：分析步总时间s
      DIMENSION COORDS(3),FLUX(2),TIME(2)
      CHARACTER*80 SNAME
C     v:热源移动速度mm/s
      v=8
C     q:热源功率mJ/s
      q=1500000
      d=v*TIME(2)

C     热源当前坐标mm
      x=COORDS(1)
      y=COORDS(2)
      z=COORDS(3)

C     初始热源中心mm
      x0=0
      y0=6
      z0=0

C     热源形状参数mm
      a=3.8
      b=4.2
      c=3.9
      PI=3.1415926

      heat=6*sqrt(3.0)*q/(a*b*c*PI*sqrt(PI))

      heat=6*sqrt(3.0)*q/(a*b*c*PI*sqrt(PI))
      shape=exp(-3*(x-x0)**2/a**2-3*(y-y0)**2/b**2-3*(z-z0-d)**2/c**2)
C     JLTYP＝1，表示为体热源
      JLTYP=1
      IF (KSTEP.eq.one) THEN
          FLUX(1)=heat*shape
      ENDIF
      RETURN
      END